32. 长三角上市公司因子聚类（信用卡样例对照）

from pathlib import Path  # 导入路径工具以定位真实因子与公司快照
import numpy as np  # 导入数值工具以对数化市值并构造随机基线
import pandas as pd  # 导入表格工具以筛选公司截面
data_root=Path('/home/ubuntu/r2_data_mount/data/stock')  # 绑定课堂股票数据目录
factor_path=data_root/'valuation_factors_quarterly_15_years.h5'  # 定位季度估值因子面板
basic_path=data_root/'stock_basic_data.h5'  # 定位公司地区与行业主表
if not factor_path.exists() or not basic_path.exists(): raise FileNotFoundError('未找到课程数据文件，请从课程数据下载入口获取并核对文件位置')  # 规定资产缺失时终止聚类
valuation_panel=pd.read_hdf(factor_path,key='valuation_factors')  # 读取真实A股季度估值因子
company_basic=pd.read_hdf(basic_path,key='stock_basic_info')  # 读取公司的省份与行业信息
company_basic['province']=company_basic['province'].astype('string')  # 将省份标签统一为字符串类型
delta_company_codes=set(company_basic.loc[company_basic['province'].isin(['上海市','浙江省','江苏省','安徽省']),'order_book_id'])  # 提取长三角公司代码
factor_columns=['ep_ratio_ttm','book_to_market_ratio_lf','dividend_yield_ttm','market_cap']  # 固定盈利、价值、股息与规模因子
quarter_list=sorted(valuation_panel.index.get_level_values('date').unique())  # 排序面板内全部报告期
quarter_now,quarter_prev=quarter_list[-1],quarter_list[-2]  # 锁定最新与上一报告期

def build_quarter_features(quarter):  # 定义长三角公司截面标准化函数
    cross_section=valuation_panel.loc[(slice(None),quarter),factor_columns].dropna().reset_index('date',drop=True)  # 提取当季完整因子截面
    cross_section=cross_section.loc[cross_section.index.isin(delta_company_codes)]  # 仅保留长三角公司
    cross_section=cross_section.assign(log_market_cap=np.log(cross_section['market_cap'])).drop(columns='market_cap')  # 对数化市值并删除原尺度列
    return (cross_section-cross_section.mean())/cross_section.std()  # 在当期截面内拟合z分数尺度
x_now=build_quarter_features(quarter_now)  # 构造最新报告期特征矩阵
x_prev=build_quarter_features(quarter_prev)  # 构造上一报告期特征矩阵
common_codes=x_now.index.intersection(x_prev.index)  # 对齐两季共同公司以核对稳定性
if len(common_codes)<20: raise ValueError('当前样本不符合模型使用条件')  # 共同公司不足时终止聚类采用

from sklearn.cluster import KMeans  # 导入K-Means以拟合候选簇数
from sklearn.metrics import silhouette_score,adjusted_rand_score  # 导入分类或聚类指标以评价样本外结果
silhouette_scores={cluster_count:silhouette_score(x_now,KMeans(cluster_count,n_init=20,random_state=42).fit_predict(x_now)) for cluster_count in range(2,7)}  # 生成silhouette_scores以承接当前分析步骤
best_cluster_count=max(silhouette_scores,key=silhouette_scores.get)  # 选择轮廓系数最高的候选簇数
current_labels=KMeans(best_cluster_count,n_init=20,random_state=42).fit_predict(x_now)  # 选择轮廓系数最高的候选簇数
previous_labels=KMeans(best_cluster_count,n_init=20,random_state=42).fit_predict(x_prev.loc[common_codes])  # 选择轮廓系数最高的候选簇数
current_common_labels=pd.Series(current_labels,index=x_now.index).loc[common_codes]  # 按共同主体对齐当前标签
stability_ari=adjusted_rand_score(previous_labels,current_common_labels); cluster_shares=pd.Series(current_labels).value_counts(normalize=True)  # 计算各簇样本占比以检查可运营规模
random_baseline=adjusted_rand_score(previous_labels,np.random.default_rng(42).permutation(current_common_labels))  # 生成固定置换基线以判断稳定性是否超越随机
adoption_conditions_met=silhouette_scores[best_cluster_count]>=.20 and stability_ari>=.60 and cluster_shares.min()>=.10 and stability_ari>random_baseline  # 选择轮廓系数最高的候选簇数
print({'K':best_cluster_count,'silhouette':silhouette_scores[best_cluster_count],'ARI':stability_ari,'shares':cluster_shares.to_dict(),'random':random_baseline,'conditions_met':adoption_conditions_met})  # 输出聚类质量和稳定性比较

target_industry='汽车制造业'  # 事先设定另一中国行业作为公司层级新案例对象
industry_codes=set(company_basic.loc[company_basic['industry_name'].eq(target_industry),'order_book_id'])  # 从真实公司主表取得行业证券代码
def build_industry_features(quarter):  # 按与主案例一致的因子口径构造行业截面
    industry_cross_section=valuation_panel.loc[(slice(None),quarter),factor_columns].dropna().reset_index('date',drop=True)  # 提取指定报告期完整因子
    industry_cross_section=industry_cross_section[industry_cross_section.index.isin(industry_codes)]  # 仅保留事先设定行业公司
    industry_cross_section=industry_cross_section.assign(log_market_cap=np.log(industry_cross_section['market_cap'])).drop(columns='market_cap')  # 采用与主案例一致的市值变换
    return (industry_cross_section-industry_cross_section.mean())/industry_cross_section.std()  # 在新案例行业内重新拟合标准化
industry_now=build_industry_features(quarter_now)  # 构造当前季度新案例行业截面
industry_prev=build_industry_features(quarter_prev)  # 构造上一季度新案例行业截面
industry_common_codes=industry_now.index.intersection(industry_prev.index)  # 对齐两季共同公司
if len(industry_common_codes)<20: raise ValueError('当前样本不符合模型使用条件')  # 样本不足时停止聚类采用
industry_labels_now=KMeans(best_cluster_count,n_init=20,random_state=42).fit_predict(industry_now)  # 用相同K拟合当前行业截面
industry_labels_prev=KMeans(best_cluster_count,n_init=20,random_state=42).fit_predict(industry_prev.loc[industry_common_codes])  # 用相同K拟合上一季度共同样本
industry_common_labels=pd.Series(industry_labels_now,index=industry_now.index).loc[industry_common_codes]  # 对齐当前季度共同公司标签
industry_ari=adjusted_rand_score(industry_labels_prev,industry_common_labels)  # 核对新案例行业跨期稳定性
industry_shares=pd.Series(industry_labels_now).value_counts(normalize=True)  # 计算新案例行业最小簇占比
industry_conditions_met=industry_ari>=.60 and industry_shares.min()>=.10  # 使用相同稳定性与规模门判断是否采用
print({'industry':target_industry,'rows_now':len(industry_now),'common_rows':len(industry_common_codes),'ARI':industry_ari,'shares':industry_shares.to_dict(),'conditions_met':industry_conditions_met})  # 输出新案例的样本量、稳定性和组别占比